path_proj = here::here()
path_source = file.path(path_proj, "source")

source(file.path(path_source, "simulation", "simulations_functions_final.R"))
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:reshape2':
## 
##     smiths
source(file.path(path_source, "functions", "plot_function.R"))
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
source(file.path(path_source, "functions", "fit_function.R"))
source(file.path(path_source, "functions", "table_function.R"))

# place for draws
# mac
posterior_draws_path = file.path(Sys.getenv("HOME"), "Desktop", "draws", "testEach")

#Windows
#posterior_draws_path = file.path(Sys.getenv("USERPROFILE"), "Desktop", "draws", "testEach")

#data path
data_save_path = file.path(path_proj, "data", "fitted_model", "simulation", "4. b_ou")
#models
q_constant <- file.path(path_proj, "source", "models", 
                      "q-constant.stan")
b_constant <- file.path(path_proj, "source", "models", 
                     "b-constant.stan")
b_rw <- file.path(path_proj, "source", "models", 
                     "b-rw1.stan")
b_ou <-  file.path(path_proj, "source", "models", 
                      "b-ou.stan")

compiled_models <- list(
  q_constant = cmdstan_model(q_constant),
  b_constant = cmdstan_model(b_constant),
  b_rw = cmdstan_model(b_rw),
  b_ou = cmdstan_model(b_ou)
)

models_to_use <- c("q_constant", "b_constant", "b_rw", "b_ou")

simulate data - FR

setting

###### setting #####
seed <- 123
set.seed(seed)

# data
alpha_increase_seq_1 <- seq(10, 750, by = 30)
alpha_decrease_seq_1 <- seq(750, 10, by = -30)
alpha_lamb =  c( rep(10,5), alpha_increase_seq_1 + rnorm(alpha_increase_seq_1,10,10), 
                 alpha_decrease_seq_1 + rnorm(alpha_decrease_seq_1,10,10),
                 rep(10,5))
beta_lamb = 0.5
T = 60
# reprot delay
D <- 15;

# Time period for checking
D_check <- 5

first_date <- as.Date("2024-01-01")

scoreRange <- c(first_date+days(9), first_date+days(19), first_date+days(29),
                first_date+days(39), first_date+days(49))
params_b_ou_FR <- list(
  data = list(
    alpha_lamb = alpha_lamb,
    beta_lamb  = beta_lamb,
    T       = T,
    date_start = as.Date("2024-01-01"),
    D = D
  ),
  q_model = list(
    method        = "b_ou",
    method_params = list(theta_logb = 0.3, mu_logb = log(0.7), init_logb = log(0.7), sigma_logb = 0.2,
                         theta_logitphi = 0.3, mu_logitphi = 1, init_logitphi = 1, sigma_logitphi = 0.2)
  )
)
b_ou_FR <- simulateData(params_b_ou_FR)
## [1] "here we go"
par(mfrow = c(2, 1))
plot(b_ou_FR$b, pch = 19, type = "b")
plot(b_ou_FR$phi, pch = 19, type = "b")

par(mfrow = c(1, 1))
matplot(t(b_ou_FR$q), type = "l", lty = 1, ylim = c(0, 1))

# ou_NFR
params_b_ou_NFR <- list(
  data = list(
    alpha_lamb = alpha_lamb,
    beta_lamb  = beta_lamb,
    T       = T,
    date_start = as.Date("2024-01-01"),
    D = D
  ),
  q_model = list(
    method        = "b_ou",
    method_params = list(theta_logb = 0.2, mu_logb = log(0.1), init_logb = log(0.1), sigma_logb = 0.15,
                         theta_logitphi = 0.2, mu_logitphi = 1.5, init_logitphi = 1.5, sigma_logitphi = 0.15)
  )
)

b_ou_NFR <- simulateData(params_b_ou_NFR)
## [1] "here we go"
par(mfrow = c(2, 1))
plot(b_ou_NFR$b, pch = 19, type = "b")
plot(b_ou_NFR$phi, pch = 19, type = "b")

par(mfrow = c(1, 1))
matplot(t(b_ou_NFR$q), type = "l", lty = 1, ylim = c(0, 1))

expplot

# exploritary analysis
page_num <- ceiling(nrow(b_ou_FR$case_reported_cumulated)/16)
exp_plot_ou <- fit_exp_plot(b_ou_FR$case_reported_cumulated,ncol = 4, nrow = 4, page = c(1:page_num), if_fit = T)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
print(exp_plot_ou)
## $plots
## $plots[[1]]

## 
## $plots[[2]]

## 
## $plots[[3]]

## 
## $plots[[4]]

## 
## 
## $coefficients
##            b          phi
## 1  1.0184186  0.253431592
## 2  0.6509694  0.546903210
## 3  0.7644895  0.136673976
## 4  1.2869640  0.203627561
## 5  0.3912118  0.300133052
## 6  0.7621946  0.007533969
## 7  0.6660028  0.171683278
## 8  0.5773877  0.254000196
## 9  0.4946644  0.268279415
## 10 0.5883842  0.239827799
## 11 0.7630322  0.202700208
## 12 0.6623841  0.304778809
## 13 0.9003512  0.321959789
## 14 0.6869329  0.280528247
## 15 0.6180896  0.314059378
## 16 0.6570050  0.316283697
## 17 0.6832462  0.268327876
## 18 0.7039898  0.258409132
## 19 0.6313040  0.252188622
## 20 0.8399478  0.240887664
## 21 0.7790189  0.198963504
## 22 0.8458926  0.198780186
## 23 0.8819860  0.236302753
## 24 1.0356185  0.272928726
## 25 1.4930619  0.213407152
## 26 1.0831836  0.271403455
## 27 0.8463732  0.250994725
## 28 0.8031341  0.260904728
## 29 0.6290154  0.284322147
## 30 0.5291156  0.333024875
## 31 0.5792335  0.284138676
## 32 0.5263670  0.253988563
## 33 0.4082156  0.272677492
## 34 0.5469907  0.268459929
## 35 0.5964300  0.272924992
## 36 0.5389308  0.311407775
## 37 0.5722262  0.347599713
## 38 0.5731940  0.333035066
## 39 0.8474409  0.362714288
## 40 0.8816275  0.302867356
## 41 0.6565743  0.316318277
## 42 0.8850903  0.283031911
## 43 0.8945391  0.265147315
## 44 0.5466306  0.233898490
## 45 0.4094268  0.220531416
## 46 0.7308413  0.306687247
## 47 0.8409546  0.305236260
## 48 0.5549634  0.334415079
## 49 0.5547306  0.357562807
## 50 0.3802698  0.281724332
## 51 0.5995839  0.325938598
## 52 0.8741929  0.291401039
## 53 0.7899340  0.326956538
## 54 0.7464737  0.338130582
## 55 0.9055233  0.343743089
## 56 1.1627079  0.426495394
## 57 1.3283720  0.441025876
## 58 0.8215360  0.055076447
## 59 0.6494544  0.294071003
## 60 0.9061664 -0.039585126
exp_plot_ou$coefficients
##            b          phi
## 1  1.0184186  0.253431592
## 2  0.6509694  0.546903210
## 3  0.7644895  0.136673976
## 4  1.2869640  0.203627561
## 5  0.3912118  0.300133052
## 6  0.7621946  0.007533969
## 7  0.6660028  0.171683278
## 8  0.5773877  0.254000196
## 9  0.4946644  0.268279415
## 10 0.5883842  0.239827799
## 11 0.7630322  0.202700208
## 12 0.6623841  0.304778809
## 13 0.9003512  0.321959789
## 14 0.6869329  0.280528247
## 15 0.6180896  0.314059378
## 16 0.6570050  0.316283697
## 17 0.6832462  0.268327876
## 18 0.7039898  0.258409132
## 19 0.6313040  0.252188622
## 20 0.8399478  0.240887664
## 21 0.7790189  0.198963504
## 22 0.8458926  0.198780186
## 23 0.8819860  0.236302753
## 24 1.0356185  0.272928726
## 25 1.4930619  0.213407152
## 26 1.0831836  0.271403455
## 27 0.8463732  0.250994725
## 28 0.8031341  0.260904728
## 29 0.6290154  0.284322147
## 30 0.5291156  0.333024875
## 31 0.5792335  0.284138676
## 32 0.5263670  0.253988563
## 33 0.4082156  0.272677492
## 34 0.5469907  0.268459929
## 35 0.5964300  0.272924992
## 36 0.5389308  0.311407775
## 37 0.5722262  0.347599713
## 38 0.5731940  0.333035066
## 39 0.8474409  0.362714288
## 40 0.8816275  0.302867356
## 41 0.6565743  0.316318277
## 42 0.8850903  0.283031911
## 43 0.8945391  0.265147315
## 44 0.5466306  0.233898490
## 45 0.4094268  0.220531416
## 46 0.7308413  0.306687247
## 47 0.8409546  0.305236260
## 48 0.5549634  0.334415079
## 49 0.5547306  0.357562807
## 50 0.3802698  0.281724332
## 51 0.5995839  0.325938598
## 52 0.8741929  0.291401039
## 53 0.7899340  0.326956538
## 54 0.7464737  0.338130582
## 55 0.9055233  0.343743089
## 56 1.1627079  0.426495394
## 57 1.3283720  0.441025876
## 58 0.8215360  0.055076447
## 59 0.6494544  0.294071003
## 60 0.9061664 -0.039585126
exp_b_data_ou<- data.frame( date = as.Date(rownames(b_ou_FR$case_reported_cumulated)),
                          b = exp_plot_ou$coefficients$b)

exp_b_plot_ou <- ggplot(exp_b_data_ou, aes(x = date, y = b)) +
  geom_point(color = "black", size = 1.5) +       
  geom_smooth(method = "loess", se = TRUE,        
              color = "blue", fill = "grey", alpha = 0.5) +
  theme_minimal() +
  labs(x = NULL, y = "Y", title = "Smoothed Curve of parameter b")

print(exp_b_plot_ou)
## `geom_smooth()` using formula = 'y ~ x'

# exploritary analysis
page_num <- ceiling(nrow(b_ou_NFR$case_reported_cumulated)/16)
exp_plot_ou <- fit_exp_plot(b_ou_NFR$case_reported_cumulated,ncol = 4, nrow = 4, page = c(1:page_num), if_fit = T)
print(exp_plot_ou)
## $plots
## $plots[[1]]

## 
## $plots[[2]]

## 
## $plots[[3]]

## 
## $plots[[4]]

## 
## 
## $coefficients
##            b         phi
## 1  0.1607346  0.08526495
## 2  0.1387379  0.39334920
## 3  0.1881441  0.30373372
## 4  0.1005666  0.23516450
## 5  0.1931866  0.05402949
## 6  0.1922913  0.48818141
## 7  0.1638292  0.17874300
## 8  0.1546258  0.22823144
## 9  0.1463539  0.25607718
## 10 0.1533422  0.24070483
## 11 0.1691421  0.23803313
## 12 0.1565306  0.19818969
## 13 0.1556989  0.22650387
## 14 0.1620785  0.18961150
## 15 0.1713553  0.21191501
## 16 0.1803743  0.14048080
## 17 0.1692417  0.17785194
## 18 0.1847041  0.10323804
## 19 0.1889704  0.09742637
## 20 0.1676988  0.11164840
## 21 0.1624350  0.09863313
## 22 0.1811429  0.08812489
## 23 0.1678583  0.13704453
## 24 0.1649921  0.11801358
## 25 0.1659259  0.11456425
## 26 0.1636334  0.15802306
## 27 0.1595850  0.14412192
## 28 0.1633322  0.19151782
## 29 0.1601417  0.19807313
## 30 0.1761062  0.21671626
## 31 0.1720761  0.17107231
## 32 0.1948724  0.20224538
## 33 0.1624139  0.16540786
## 34 0.1649930  0.14688458
## 35 0.1687533  0.19983575
## 36 0.1853506  0.14074984
## 37 0.1802496  0.20171508
## 38 0.1746101  0.21637612
## 39 0.1834760  0.16829431
## 40 0.1660187  0.18218919
## 41 0.1955013  0.16014742
## 42 0.1834835  0.12953189
## 43 0.2307871  0.20862896
## 44 0.1635083  0.16687631
## 45 0.1733182  0.16030699
## 46 0.1721403  0.17139965
## 47 0.1729932  0.12502694
## 48 0.1607783  0.18134152
## 49 0.1737043  0.13383785
## 50 0.1866806  0.15535871
## 51 0.1770962  0.13958139
## 52 0.1622209  0.19588355
## 53 0.1645726  0.08569237
## 54 0.1476060  0.22149815
## 55 0.1858383  0.22147106
## 56 0.1692664  0.21831461
## 57 0.2532632  0.27274666
## 58 0.2076270  0.15474529
## 59 0.1582585 -0.10456069
## 60 0.1893347  0.03540411
exp_plot_ou$coefficients
##            b         phi
## 1  0.1607346  0.08526495
## 2  0.1387379  0.39334920
## 3  0.1881441  0.30373372
## 4  0.1005666  0.23516450
## 5  0.1931866  0.05402949
## 6  0.1922913  0.48818141
## 7  0.1638292  0.17874300
## 8  0.1546258  0.22823144
## 9  0.1463539  0.25607718
## 10 0.1533422  0.24070483
## 11 0.1691421  0.23803313
## 12 0.1565306  0.19818969
## 13 0.1556989  0.22650387
## 14 0.1620785  0.18961150
## 15 0.1713553  0.21191501
## 16 0.1803743  0.14048080
## 17 0.1692417  0.17785194
## 18 0.1847041  0.10323804
## 19 0.1889704  0.09742637
## 20 0.1676988  0.11164840
## 21 0.1624350  0.09863313
## 22 0.1811429  0.08812489
## 23 0.1678583  0.13704453
## 24 0.1649921  0.11801358
## 25 0.1659259  0.11456425
## 26 0.1636334  0.15802306
## 27 0.1595850  0.14412192
## 28 0.1633322  0.19151782
## 29 0.1601417  0.19807313
## 30 0.1761062  0.21671626
## 31 0.1720761  0.17107231
## 32 0.1948724  0.20224538
## 33 0.1624139  0.16540786
## 34 0.1649930  0.14688458
## 35 0.1687533  0.19983575
## 36 0.1853506  0.14074984
## 37 0.1802496  0.20171508
## 38 0.1746101  0.21637612
## 39 0.1834760  0.16829431
## 40 0.1660187  0.18218919
## 41 0.1955013  0.16014742
## 42 0.1834835  0.12953189
## 43 0.2307871  0.20862896
## 44 0.1635083  0.16687631
## 45 0.1733182  0.16030699
## 46 0.1721403  0.17139965
## 47 0.1729932  0.12502694
## 48 0.1607783  0.18134152
## 49 0.1737043  0.13383785
## 50 0.1866806  0.15535871
## 51 0.1770962  0.13958139
## 52 0.1622209  0.19588355
## 53 0.1645726  0.08569237
## 54 0.1476060  0.22149815
## 55 0.1858383  0.22147106
## 56 0.1692664  0.21831461
## 57 0.2532632  0.27274666
## 58 0.2076270  0.15474529
## 59 0.1582585 -0.10456069
## 60 0.1893347  0.03540411
exp_b_data_ou<- data.frame( date = as.Date(rownames(b_ou_NFR$case_reported_cumulated)),
                          b = exp_plot_ou$coefficients$b)

exp_b_plot_ou <- ggplot(exp_b_data_ou, aes(x = date, y = b)) +
  geom_point(color = "black", size = 1.5) +       
  geom_smooth(method = "loess", se = TRUE,        
              color = "blue", fill = "grey", alpha = 0.5) +
  theme_minimal() +
  labs(x = NULL, y = "Y", title = "Smoothed Curve of parameter b")

print(exp_b_plot_ou)
## `geom_smooth()` using formula = 'y ~ x'

fit model

out_b_ou_FR <-
  nowcasting_moving_window(b_ou_FR$case_reported_cumulated, scoreRange = scoreRange,
                          case_true = b_ou_FR$case_true,
                          start_date = as.Date("2024-01-01"),
                          D = D, seeds = seed,
                          methods =models_to_use,
                          compiled_models = compiled_models,
                          iter_sampling = 2000, iter_warmup = 1000, refresh = 0,
                          num_chains = 3, suppress_output = T,
                          posterior_draws_path = file.path(posterior_draws_path, "b_ou")
                          )
## ====================
## now=2024-01-10 (1/5)
## ====================
## Warning in nowcasting_moving_window(b_ou_FR$case_reported_cumulated, scoreRange
## = scoreRange, : The number of rows of the input data is smaller than number of
## max delay D, which might cause inaccuracy.
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 0.7 seconds.
## Chain 2 finished in 0.6 seconds.
## Chain 3 finished in 0.6 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 0.6 seconds.
## Total execution time: 2.3 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 0.2 seconds.
## Chain 2 finished in 0.2 seconds.
## Chain 3 finished in 0.2 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 1.1 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 2.4 seconds.
## Chain 2 finished in 2.7 seconds.
## Chain 3 finished in 1.8 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 2.3 seconds.
## Total execution time: 7.2 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 3.2 seconds.
## Chain 2 finished in 3.0 seconds.
## Chain 3 finished in 1.3 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 2.5 seconds.
## Total execution time: 7.8 seconds.
## 
## ====================
## now=2024-01-20 (2/5)
## ====================
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 1.6 seconds.
## Chain 2 finished in 1.4 seconds.
## Chain 3 finished in 1.5 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 1.5 seconds.
## Total execution time: 4.9 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 0.4 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 3 finished in 0.4 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 1.4 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 3.5 seconds.
## Chain 2 finished in 3.3 seconds.
## Chain 3 finished in 3.5 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 3.4 seconds.
## Total execution time: 10.6 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 3.3 seconds.
## Chain 2 finished in 5.4 seconds.
## Chain 3 finished in 3.1 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 3.9 seconds.
## Total execution time: 12.1 seconds.
## 
## ====================
## now=2024-01-30 (3/5)
## ====================
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 2.2 seconds.
## Chain 2 finished in 2.0 seconds.
## Chain 3 finished in 1.9 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 2.0 seconds.
## Total execution time: 6.4 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 0.6 seconds.
## Chain 2 finished in 0.6 seconds.
## Chain 3 finished in 0.6 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 0.6 seconds.
## Total execution time: 2.0 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 4.3 seconds.
## Chain 2 finished in 4.3 seconds.
## Chain 3 finished in 4.5 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 4.4 seconds.
## Total execution time: 13.4 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 3.5 seconds.
## Chain 2 finished in 3.9 seconds.
## Chain 3 finished in 3.3 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 3.6 seconds.
## Total execution time: 11.1 seconds.
## 
## ====================
## now=2024-02-09 (4/5)
## ====================
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 2.6 seconds.
## Chain 2 finished in 2.7 seconds.
## Chain 3 finished in 2.8 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 2.7 seconds.
## Total execution time: 8.5 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 0.8 seconds.
## Chain 2 finished in 0.8 seconds.
## Chain 3 finished in 0.7 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 0.8 seconds.
## Total execution time: 2.6 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 5.2 seconds.
## Chain 2 finished in 6.4 seconds.
## Chain 3 finished in 6.7 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 6.1 seconds.
## Total execution time: 18.5 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 5.3 seconds.
## Chain 2 finished in 6.0 seconds.
## Chain 3 finished in 5.3 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 5.6 seconds.
## Total execution time: 17.0 seconds.
## 
## ====================
## now=2024-02-19 (5/5)
## ====================
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 3.6 seconds.
## Chain 2 finished in 3.5 seconds.
## Chain 3 finished in 3.8 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 3.6 seconds.
## Total execution time: 11.3 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 1.1 seconds.
## Chain 2 finished in 1.0 seconds.
## Chain 3 finished in 1.0 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 1.0 seconds.
## Total execution time: 3.4 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 5.9 seconds.
## Chain 2 finished in 7.7 seconds.
## Chain 3 finished in 6.1 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 6.6 seconds.
## Total execution time: 20.0 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 6.9 seconds.
## Chain 2 finished in 5.8 seconds.
## Chain 3 finished in 6.0 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 6.2 seconds.
## Total execution time: 19.2 seconds.
save(out_b_ou_FR, file = file.path(data_save_path, "FR_b_ou.RData"))
#load(file.path(data_save_path, "FR_b_ou.RData"))

results_b_ou_FR <- nowcasts_table(out_b_ou_FR, D = D, report_unit = "day", 
                          methods = models_to_use)

plots_b_ou_FR <- nowcasts_plot(results_b_ou_FR, D = D, report_unit = "day", methods = models_to_use)
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
for (i in 1:length(results_b_ou_FR)) {
  print(calculate_metrics(results_b_ou_FR[[i]],
                          methods = models_to_use))
}
##    RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 62.94 59.77 43.77 59.08         106.61           0.7 q_constant
## 2 15.73  9.25  9.35  8.04          50.00           1.0 b_constant
## 3 16.49  8.63  8.98  7.17          62.40           1.0       b_rw
## 4 17.59  9.89 10.22  8.49          63.10           1.0       b_ou
##    RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 53.88 12.68 39.27 12.00         104.11          0.65 q_constant
## 2 29.59  4.02 11.49  2.55          77.56          1.00 b_constant
## 3 31.99  4.19 13.48  2.64          96.11          1.00       b_rw
## 4 26.08  3.40  9.76  2.00          93.46          1.00       b_ou
##     RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 110.57  8.21 47.47 4.94         110.07          0.87 q_constant
## 2  87.04  6.23 28.23 2.95         103.51          0.90 b_constant
## 3  36.53  2.74 14.45 1.53         135.61          1.00       b_rw
## 4  51.51  3.61 14.88 1.52         135.18          1.00       b_ou
##    RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 48.02  4.24 24.17 2.48         112.53          0.95 q_constant
## 2 48.38  4.18 24.53 2.54         110.31          0.95 b_constant
## 3 29.25  2.78  9.62 1.18         138.78          1.00       b_rw
## 4 31.83  3.00  9.20 1.14         139.44          1.00       b_ou
##    RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 26.15  3.47 17.15 2.40         108.01          0.98 q_constant
## 2 25.02  3.28 15.20 2.10         106.35          0.98 b_constant
## 3 17.36  3.33  7.72 1.51         115.93          0.98       b_rw
## 4 14.54  2.68  6.42 1.26         114.83          1.00       b_ou
for (i in 1:length(results_b_ou_FR)) {
  print(calculate_metrics(data.table::last(results_b_ou_FR[[i]],D_check),
                          methods = models_to_use))
}
##    RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 88.34 53.64 77.03 53.49         185.02           0.4 q_constant
## 2 22.22 11.23 17.83 10.96          82.61           1.0 b_constant
## 3 23.30 10.74 17.21  9.87         108.01           1.0       b_rw
## 4 24.84 12.02 19.49 11.47         109.20           1.0       b_ou
##    RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 80.13 11.96 72.25 10.52         217.63           0.4 q_constant
## 2 58.33  6.90 36.94  4.46         160.20           1.0 b_constant
## 3 63.35  7.51 45.26  5.52         231.62           1.0       b_rw
## 4 51.73  6.11 32.31  3.85         220.82           1.0       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 255.64 17.62 175.11 12.47         244.02           0.6 q_constant
## 2 208.02 14.21 121.62  8.53         228.02           0.6 b_constant
## 3  88.70  6.15  72.48  5.13         416.23           1.0       b_rw
## 4 125.84  8.56  76.90  5.32         413.24           1.0       b_ou
##     RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 106.92  9.85 75.00 6.84         197.61           0.8 q_constant
## 2 101.09  9.28 74.16 6.74         190.42           0.8 b_constant
## 3  79.97  7.37 54.59 4.98         385.43           1.0       b_rw
## 4  87.98  8.08 52.54 4.76         389.63           1.0       b_ou
##    RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 26.02  5.11 22.22 4.39         126.60             1 q_constant
## 2 21.26  4.19 16.71 3.28         124.83             1 b_constant
## 3 42.47  8.93 36.21 7.40         199.22             1       b_rw
## 4 32.44  6.66 28.06 5.67         187.41             1       b_ou
print(plots_b_ou_FR)
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

out_b_ou_NFR <-
  nowcasting_moving_window(b_ou_NFR$case_reported_cumulated, scoreRange = scoreRange,
                          case_true = b_ou_NFR$case_true,
                          start_date = as.Date("2024-01-01"),
                          D = D, seeds = seed,
                          methods =models_to_use,
                          compiled_models = compiled_models,
                          iter_sampling = 2000, iter_warmup = 1000, refresh = 0,
                          num_chains = 3, suppress_output = T,
                          posterior_draws_path = file.path(posterior_draws_path, "b_ou")
                          )
## ====================
## now=2024-01-10 (1/5)
## ====================
## Warning in nowcasting_moving_window(b_ou_NFR$case_reported_cumulated,
## scoreRange = scoreRange, : The number of rows of the input data is smaller than
## number of max delay D, which might cause inaccuracy.
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 0.7 seconds.
## Chain 2 finished in 0.7 seconds.
## Chain 3 finished in 0.7 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 0.7 seconds.
## Total execution time: 2.4 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 0.3 seconds.
## Chain 2 finished in 0.2 seconds.
## Chain 3 finished in 0.3 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 0.3 seconds.
## Total execution time: 1.1 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 3.2 seconds.
## Chain 2 finished in 2.1 seconds.
## Chain 3 finished in 2.7 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 2.7 seconds.
## Total execution time: 8.3 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 2.2 seconds.
## Chain 2 finished in 3.9 seconds.
## Chain 3 finished in 2.6 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 2.9 seconds.
## Total execution time: 9.1 seconds.
## 
## ====================
## now=2024-01-20 (2/5)
## ====================
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 1.4 seconds.
## Chain 2 finished in 1.5 seconds.
## Chain 3 finished in 1.4 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 1.4 seconds.
## Total execution time: 4.7 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 0.9 seconds.
## Chain 2 finished in 0.8 seconds.
## Chain 3 finished in 0.9 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 0.9 seconds.
## Total execution time: 2.9 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 5.6 seconds.
## Chain 2 finished in 5.2 seconds.
## Chain 3 finished in 6.1 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 5.6 seconds.
## Total execution time: 17.2 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 3.9 seconds.
## Chain 2 finished in 4.8 seconds.
## Chain 3 finished in 5.2 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 4.6 seconds.
## Total execution time: 14.1 seconds.
## 
## ====================
## now=2024-01-30 (3/5)
## ====================
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 1.6 seconds.
## Chain 2 finished in 1.7 seconds.
## Chain 3 finished in 1.7 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 1.7 seconds.
## Total execution time: 5.4 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 1.1 seconds.
## Chain 2 finished in 1.1 seconds.
## Chain 3 finished in 1.1 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 1.1 seconds.
## Total execution time: 3.6 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 9.3 seconds.
## Chain 2 finished in 8.9 seconds.
## Chain 3 finished in 11.1 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 9.8 seconds.
## Total execution time: 29.6 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 6.5 seconds.
## Chain 2 finished in 8.1 seconds.
## Chain 3 finished in 6.8 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 7.1 seconds.
## Total execution time: 21.8 seconds.
## 
## ====================
## now=2024-02-09 (4/5)
## ====================
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 2.1 seconds.
## Chain 2 finished in 1.9 seconds.
## Chain 3 finished in 1.7 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 1.9 seconds.
## Total execution time: 6.1 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 1.7 seconds.
## Chain 2 finished in 1.5 seconds.
## Chain 3 finished in 1.8 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 1.7 seconds.
## Total execution time: 5.4 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 13.5 seconds.
## Chain 2 finished in 12.9 seconds.
## Chain 3 finished in 14.4 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 13.6 seconds.
## Total execution time: 41.2 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 12.7 seconds.
## Chain 2 finished in 10.9 seconds.
## Chain 3 finished in 11.3 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 11.6 seconds.
## Total execution time: 35.3 seconds.
## 
## ====================
## now=2024-02-19 (5/5)
## ====================
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 2.1 seconds.
## Chain 2 finished in 2.5 seconds.
## Chain 3 finished in 2.7 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 2.4 seconds.
## Total execution time: 7.6 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 1.6 seconds.
## Chain 2 finished in 1.6 seconds.
## Chain 3 finished in 1.7 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 1.6 seconds.
## Total execution time: 5.2 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 13.3 seconds.
## Chain 2 finished in 15.8 seconds.
## Chain 3 finished in 17.7 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 15.6 seconds.
## Total execution time: 47.1 seconds.
## 
## Running MCMC with 3 sequential chains...
## 
## Chain 1 finished in 14.9 seconds.
## Chain 2 finished in 14.0 seconds.
## Chain 3 finished in 11.5 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 13.5 seconds.
## Total execution time: 40.7 seconds.
save(out_b_ou_NFR, file = file.path(data_save_path, "NFR_b_ou.RData"))
# load(file.path(data_save_path, "NFR_b_ou.RData"))

results_b_ou_NFR <- nowcasts_table(out_b_ou_NFR, D = D, report_unit = "day",
                          methods = models_to_use)

plots_b_ou_NFR <- nowcasts_plot(results_b_ou_NFR, D = D, report_unit = "day", methods = models_to_use)

for (i in 1:length(results_b_ou_NFR)) {
  print(calculate_metrics(results_b_ou_NFR[[i]],
                          methods = models_to_use))
}
##    RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 18.24 21.00 13.27 18.51          82.00           1.0 q_constant
## 2 67.66 49.19 45.39 47.61          41.30           0.5 b_constant
## 3 70.83 50.09 47.19 48.45          46.30           0.5       b_rw
## 4 69.75 50.11 46.65 48.65          44.81           0.4       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 156.22 27.56 104.51 25.22          90.06          0.35 q_constant
## 2 159.02 27.83 106.85 25.26          88.65          0.35 b_constant
## 3 161.04 24.68 102.16 20.83         152.01          0.60       b_rw
## 4 156.40 25.17 101.42 21.88         144.82          0.65       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 162.84 22.83 125.23 21.13         103.78          0.23 q_constant
## 2 109.16 13.61  66.98 11.01         120.85          0.57 b_constant
## 3  80.59 11.18  53.03  9.20         254.67          0.90       b_rw
## 4  71.68 11.60  48.32  9.22         209.59          0.87       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 151.22 18.90 118.55 16.84         109.21          0.25 q_constant
## 2 102.75 10.97  67.18  9.18         132.96          0.65 b_constant
## 3  68.73  9.38  46.54  7.33         233.12          0.92       b_rw
## 4  57.76  9.24  41.09  7.06         223.87          0.95       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 153.79 19.81 127.67 18.25         101.63          0.16 q_constant
## 2  77.65 10.64  54.51  8.37         115.57          0.70 b_constant
## 3  62.07 11.23  45.61  8.27         171.99          0.90       b_rw
## 4  60.71 10.59  41.41  7.62         174.75          0.88       b_ou
for (i in 1:length(results_b_ou_NFR)) {
  print(calculate_metrics(data.table::last(results_b_ou_NFR[[i]],D_check),
                          methods = models_to_use))
}
##    RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 25.59 22.50 23.80 20.53         140.80           1.0 q_constant
## 2 95.29 51.40 82.63 50.01          68.20           0.2 b_constant
## 3 99.79 53.35 86.27 52.00          77.80           0.2       b_rw
## 4 98.26 52.82 85.08 51.64          75.41           0.2       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 291.15 36.46 274.34 35.18         180.82           0.0 q_constant
## 2 295.47 37.03 278.81 35.78         177.61           0.0 b_constant
## 3 307.80 38.99 297.00 38.35         344.63           0.0       b_rw
## 4 296.19 37.33 283.66 36.52         328.44           0.2       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 196.89 15.74 184.28 14.41         247.01           0.2 q_constant
## 2 180.31 12.68 128.24  9.46         292.42           0.4 b_constant
## 3 126.37  9.56 108.97  8.35         846.54           1.0       b_rw
## 4 111.42  8.33  91.85  7.01         614.27           1.0       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1  97.43  8.65  80.56  7.26         209.21           0.6 q_constant
## 2 183.11 16.11 164.75 14.67         257.61           0.6 b_constant
## 3 103.42  9.91  82.01  7.72         625.44           1.0       b_rw
## 4  72.65  6.54  65.51  5.96         548.24           1.0       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 113.30 23.96 109.87 23.12         122.61           0.2 q_constant
## 2  51.78 11.16  42.99  9.18         143.01           0.6 b_constant
## 3 105.47 22.88 101.57 21.75         259.82           0.8       b_rw
## 4  68.27 15.08  59.52 13.09         238.20           0.8       b_ou
print(plots_b_ou_NFR)
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]